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1 Introduction 

Current density reconstructions (CDRs) of cardiac 
activation based on non-invasively obtained MCG 
(magnetocardiogram) and ECG (electrocardiogram) 
data are a new and promising tool and might support 
diagnoses in cardiology [1]. A typical result of a 
CDR is a color-coded activation map representing 
the magnitude of the current density in a volume or 
on a surface or a current dipole distribution in 3-D 
space. For statistical data analyses (e.g. in group 
studies) a method is needed which enables a 
comparison of current density distributions for 
different time points and within groups of patients 
or volunteers. This goal can be achieved by the 
proper parameterization of current density 
distributions and the application of statistical 
analyses to the parameters extracted. In this paper, 
we present a technique called as an equivalent 
ellipsoid which has been introduced previously for 
2-D planes [2]. The usefulness of our new technique 
is demonstrated through phantom data and patient 
data. 

2 Methods 

2.1 Equivalent ellipsoid 

An equivalent ellipsoid has been defined as a 3-D 
ellipsoidal object fitted to a current density 
distribution region in which the magnitude of the 
currents is above a certain threshold. Reasons for 
choosing the equivalent ellipsoid technique have 
been its geometrical simplicity (easy visualization) 
and a straightforward interpretation of lengths and 
directions of axes in the current density 
distributions. The equivalent ellipsoid is defined by 
three orthogonal semi-axes (< a=al a , b=bl b , c=cl c ), 
where 1^ denotes the unit vector connected to the 
axis £e{a,b,c}. The fitting procedure has been 
realized in following steps: 

1. The threshold T h used for marking the most 

important region in CDR has been defined as 

T h = 10 0%* ( Q max +Gnean ) /2 Gnax, where Qmax 


is the maximum and Q mean is the mean value of 
the dipole moments in the current distribution. 

2. The center of gravity of the marked region 
(COG, Pcog) has been calculated. The position 
of the COG has been used as a new origin of 
the local coordinate system for the equivalent 
ellipsoid estimation. 

3. The direction of main axis l a of the equivalent 
ellipsoid has been computed on the basis of 
three different approaches. 

(3a) First approach , \ a has been estimated on 
the basis of the normalized position of the 
dipole with the weighted longest distance from 
the COG (LD). The weighting of the distance 
has been done by the dipole moment in order to 
avoid that small and far off dipoles have a large 
influence on the definition of the main axis. 

(3b) Second approach , the normalized position 
of the CDR maximum has been applied (PM). 
(3c) Third approach uses the dominant 
direction (DD) of the region (defined by DD = 
EMi/A t , where Mi denotes the moment of i -th 
current dipole and N T is the number of dipoles 
in the thresholded region). 

4. The local coordinate system has been rotated so 
that the vector l a represents the z-axis of the 
new coordinate system. In the rotated 
coordinate system, the normalized position of 
the dipole with the maximal distance from Pcog 
on the new x-y plane has been used as the 
direction of the second semi-axis of the 
equivalent ellipsoid (1^). 

5. The direction of the third axis has been 
determined by l c =l a xl b . 

6. The average length of all axes of the equivalent 
ellipsoid in the new orthogonal coordinate 
system (l fl , l b , l c ) has been calculated itera¬ 
tively according to 

«=£Efc".l- wiihe=E M .. <d 

C? i =i i =i 

£ e (a, b, c). 



The factor g is equal to 1 and iteratively 
increased by 1% until the inside ratio (Ir) is 
above 70%. 

The inside ratio of the equivalent ellipsoid has been 
introduced in order to estimate how many current 
dipoles are located inside the equivalent ellipsoid. 
Each dipole has been associated with a voxel 
representing a small volume around the dipole. We 
have defined a voxel as an elementary cube centered 
around a current dipole with the edge length equal 
to the minimum inter distance between the current 
dipoles in the distribution. The inside ratio Ir has 
been defined as the ratio of volume of dipole voxels 
located inside the ellipsoid to the volume of all 
voxels in the thresholded distribution. 



COG 


Figure 1: a) Current density distribution; b) Thresholded 
distribution and vector l a for different approaches; c) 
Vector l b for the longest distance approach (LD); d) 
Equivalent ellipsoid estimated on the basis of the longest 
distance approach. 

In order to assess the quality of the equivalent 
ellipsoid we have introduced a goodness factor G 0 , 
which has been defined as the volume ratio of the 
dipole voxels located inside the ellipsoid and the 
equivalent ellipsoid. The equivalent ellipsoid 
technique has been tested on two problems 
presented in the following sections. 

2.2 Extended current source model 

We have used a realistic torso phantom which was 
built as a tank-like model of epoxy resin with a 
hollow space comparable to the body surface of a 
young male subject (Fig. 2a). The phantom has 
been filled with 0.9% NaCl solution which 


corresponds to a conductivity of 1.43 S/m. Inside the 
tank, an extended source model ( butterfly source) 
has been installed (Fig.2b). 


Figure 2: a) Torso phantom model; b) Butterfly extended 
source. 

The measurements have been taken in a 
magnetically shielded room (AK3b, Vacuum- 
schmelze, Hanau, Germany) at the Biomagnetic 
Center in Jena, Germany. The magnetic field has 
been recorded with a twin dewar biomagnetometer 
system (2x31 channels) with first order axial 
gradiometers (Philips, Hamburg, Germany). The 
butterfly source has been constructed of platinum 
with an extent of 5x6 cm. A sinusoidal current with 
a frequency of 25 Hz and a peak-to-peak amplitude 
of 500 juA has been impressed. 

Current density distributions have been calculated 
with the help of the CURRY software (Neuroscan, 
Inc., Sterling, VA, USA). A boundary element 
method (BEM) model has been created from 
magnetic resonance imaging (MRI) data (8796 
triangular elements) and has served as a volume 
conductor model [3]. We have used the following 
source reconstruction methods [4] (abbreviations 
used below are given in parenthesis): minimum 
norm least squares (L2), minimum L\ norm (LI), 
low resolution electromagnetic tomography (LORI, 
LOR2, where LORI denotes a L\ norm and LOR2 a 
L 2 norm applied to LORETA). The following basic 
equation has been used in all reconstructions: 

A = |M- Ljp +A 2 |Cjf" , (2) 

data term model term 

with A as the residual variance which has to be 
minimized. M is the measurement data matrix, L the 
lead field matrix, j the source vector, A the 
regularization factor, C is the weighting matrix, and 
p d and p m have been used for the norms of the data 










































and the modeling term. For the L\ norm holds 
Pd-Pm- 1, and for the L 2 norm holds p d -p m - 2. The 
current density distributions have been 
reconstructed on a multi plane 3-D grid parallel to 
the butterfly source plane (grid: 21x21x11 dipoles, 
grid spacing: d - 10 mm). 

2.3. Patient data 

Using the 62 channel biomagnetometer system 
described above, we have measured the 
magnetocardiagram of a patient aged 70 years who 
had non-sustained ventricular tachycardia that 
developed after anterior left ventricular myocardial 
infarction and apical aneurysm at the Biomagnetic 
Center Jena. We have recorded 600 s of signals at a 
sampling rate of 1000 Hz. A noise level of 50 jT has 
been estimated in both dewars. The last 40 ms of the 
bi-directional 30 Hz highpass filtered depolarization 
signal (late potentials, LP) have been used for 
inverse computations. 

A 3-D MRI data set of the chest of the patient has 
been obtained. A BEM model consisting of the left 
and right lungs as well as the outer torso surface has 
been applied for the magnetic field computations 
(forward model). We have used a triangular mesh 
with 2990 nodes and a linear potential approach for 
each triangle. The conductivity ratio of torso to 
lungs was 5 to 1. The surface of the left ventricle 
(Zv) has been segmented from the MRI data set and 
subsequently used for the restriction of the source 
space. 

Two source configurations have been investigated. 
The first configuration has consisted of dipoles 
distributed on the surface of the left ventricle (1022 
dipoles with an average dipole spacing of 4.7 mm). 
In the second one, a regular grid of dipolar sources 
(5 mm grid spacing in x, y, and z) has been placed 
into the Zv resulting in 1715 current dipoles. The 
source parameters (dipole strength and orientation) 
have been determined through a minimum norm 
least squares algorithm ( L 2 norm) for all time points. 
The equivalent ellipsoids have been estimated on 
the basis of the distribution containing the dipoles 
with the maximum strength over the last 40 ms time 
interval (LP). 

3 Results 

3.1 Extended source in torso phantom 

In rows from up to down, Fig. 3 shows both the 
results of LI, L2, LORI, and LOR2 reconstructions 


as well as the equivalent ellipsoids estimated for 
different approaches. The most focused solution has 
been obtained with the L\ norm. As expected, the 
L 2 -norm solution has been the most distributed one. 
The reconstructions calculated with LORI and 
LOR2 could have been hardly distinguished. Only 
in the case of the L\ norm the three different 
approaches (LD, DD, PM) have yielded 
significantly different ellipsoids. 



Figure 3: Results of minimum L 2 , L 1 norms and LORETA 
reconstructions on 3-D grid 21x21x11 and estimated 
equivalent ellipsoids using longest distance (LD), 
domonant direction (DD) and position of maximum (PM) 
approaches. 























































3.2 Patient data 


4 Discussion and conclusions 


Fig. 4 shows the equivalent ellipsoids found for the 
CDRs (minimum L 2 norm) on the surface of the Iv 
and on the regular 3-D grid of points located inside 
the left ventricle. 



Figure 4: Equivalent ellipsoids fitted to the reconstructed 
CDRs located on the surface of the left ventricle (left) and 
on the regular 3-D grid points inside the left ventricle 
(right). Diaphragmal view (up), apical view (middle), and 
magnified apical view for DD approach (down). 

The PM approach has yielded larger ellipsoids for 
the surface-based reconstructions than the other two 
approaches (28.4 ml versus 17.9 and 12.0 ml). LD 
and DD yield similar semi-axes lengths, while PM 
clearly differs (23.5, 17.8, 10.2), (24.8, 15.6, 7.4) 
versus (18.7, 21.2, 17.1) mm. The PM approach also 
exhibits the lowest goodness factor for the surface- 
based reconstruction (10.9%), while the DD 
approach has the best goodness factor (25.7%). The 
difference between the surface and 3-D grid-based 
COG for the thresholded region has been equal to 
4.3 mm. 


In this paper we have presented a new technique 
which enables the postprocessing of 3-D CDRs. The 
technique has been found to be useful for the 
analysis of phantom and patient data. Three 
different approaches of the algorithm have been 
tested. For 3-D grid-based CDRs all approaches 
yielded similar results. However, for a surface- 
based CDR the DD approach has been found to be 
the most appropriate one. The center of the butterfly 
source has been well estimated by the LI and L2 
methods, while both LORI and LOR2 have 
localized the center too deep. In contrary, the LORI 
and LOR2 have given a better estimation of the 
butterfly extent based on the circumscribing 
rectangle. 

One major advantage of the equivalent ellipsoid 
technique proposed is that it extracts parameters 
from CDRs which can be used for statistical 
analyses (volume of the equivalent ellipsoid, length 
and direction of the semi-axes, COG). This allows 
both comparisons of different algorithms as 
demonstrated in this paper and comparisons of 
groups of patients or volunteers. Furthermore, the 
equivalent ellipsoid provides an easy and 
straightforward visualization of the foci of even 
very complex 3-D current density distributions. 
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